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ABSTRACT 



Astrophysical fluids under the influence of magnetic fields are often subjected to 
single-fluid or two-fluid approximations. In the case of weakly ionized plasmas however, 
this can be inappropriate due to distinct responses from the multiple constituent 
species to both collisional and non-collisional forces. As a result, in dense molecular 
clouds and proto-stcllar accretion discs for instance, the conductivity of the plasma 
may be highly anisotropic leading to phenomena such as Hall and ambipolar diffusion 
strongly influencing the dynamics. 

Diffusive processes are known to restrict the stability of conventional numerical 
schemes which are not implicit in nature. Furthermore, recent work establishes that 
a large Hall term can impose an additional severe stability limit on standard explicit 
schemes. Following a previous paper which presented the one-dimensional case, we 
describe a fully three-dimensional method which relaxes the normal restrictions on 
explicit schemes for multifluid processes. This is achieved by applying the little known 
Super TimeStepping technique to the symmetric (ambipolar) component of the evo- 
lution operator for the magnetic field in the local plasma rest-frame, and the new Hall 
Diffusion Scheme to the skew-symmetric (Hall) component. 

Key words: magnetic fields - MHD - waves - methods:numerical - ISM:clouds - dust, 
extinction 



1 INTRODUCTION 

Numerical schemes used in simulations of astrophysical plas- 
mas are frequently derived from single-fluid magnetohy- 
drodynamic (MHD) models . The most common example 
of this is ideal MHD, with assumptions including infinite 
conductivity and negligible Hall current. Extended models 
within the single-fluid framework are commonly used for fi- 
nite scalar conductivity and Hall current. Furthermore, two- 
fluid models are used when the drift of a neutral component 
through the bulk plasma is considered important. With ref- 
erence to the generalized Ohm's law, we now briefly survey 
the physical motivations for departing from models based 
on ideal MHD. The discussion makes a progression through 
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We associate the multiplicity of the fluids described by a model 
to the number of fluids treated distinctly in the derived numerical 
scheme. 



various models arriving at the argument for a fully multifluid 
numerical approach to weakly ionized plasmas. 

The generalized Ohm's law for collisional gases de- 
scribes the dependencies of electric currents on the relative 
drift of charged particles due to effects both mediated by, 
and independent of, magnetic fields. In the latter case, for 
example, electron pressure can cause electrons in a local con- 
densation of gas to diffuse more quickly than ions due to 
greater thermal velocities. The resulting separation of charge 
creates an electric force coupling the ion and electron gases 
in a process known in plasma physics as ambipolar diffusion 
(Cowling 1956). In the following however, electron pressure 
is neglected under the assumptions that L c/w pe and 
L^> r c where L is the scale length of the plasma, io pc is the 
electron plasma frequency, and r c is the electron gyroradius. 
The term ambipolar diffusion is now used without ambigu- 
ity to describe an entirely different, magnetically mediated 
phenomenon of neutral drift, as more commonly discussed in 
astrophysical contexts (Mestel & Spitzer 1956; Spitzer 1978; 
and more recently, Wardle & Ng 1999). 

Defining E' as the electric field in the local rest frame 
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of the bulk plasma, and considering only effects dependent 
on the presence of a magnetic field, the generalized Ohm's 
law can be written as 

E' = a 1 J 

= roJ\\ + r u J± x B + r A J±- (1) 

In this equation er is the tensor conductivity of the plasma, 
and ro, rn, r\ are the corresponding Ohmic (field-parallel), 
Hall, and ambipolar (Pedersen) resistivities respectively. 
The explicit form of the conductivity for a weakly ionized 
plasma will be discussed in section 2, however, it is worth 
pointing out some general properties of equation (1) before 
proceeding. 

While collisions may produce rich and complex physics 
via their influence on currents, Hall diffusion can operate 
independently of collisional forces. (Note that we refer to 
the Hall term as diffusive in the sense that it contributes to 
the violation of field freezing, however, it is dispersive in na- 
ture and twists, rather than diffuses, the magnetic field.) 
Considering firstly the special case of fully ionized gases 
where r\ = ro = r rcs , the Ohmic and ambipolar terms in 
equation (f) may be combined into a single resistive term 
r les J. For L<c/u} p i, where a> p i is the ion cyclotron fre- 
quency, the greater inertia of the ions causes them to de- 
couple from the electrons (even when collisions are unim- 
portant and r rcs — > 0) and the Hall term rnJ± x B in 
equation (1) becomes significant. This regime is frequently 
approximated via the single-fluid Hall-MHD model (see, for 
example, Huba 2005; Mininni, Gomez & Mahajan 2005). 

Furthermore, when collisions are important, disparate 
resistive effects impede the flows of currents in senses both 
parallel and perpendicular to the magnetic field. In fully 
ionized plasmas, or weakly ionized plasmas where magnetic 
forces on the charged species are dominated by collisional 
drag on the neutrals, the electron drift with respect to the 
bulk plasma is fully determined by the electric current and 
a single-fluid model is tenable. Moreover, if the Hall ef- 
fect is negligible (L 3> c/Wpi), so-called resistive MHD is 
retrieved with a scalar conductivity a Tes = and corre- 
sponding Ohm's law E' = r rcs J. 

In incompletely ionized plasmas when magnetic forces 
on the charged species dominate collisional drag, ambipolar 
diffusion occurs as the charged particles remain tightly cou- 
pled to the magnetic field while drifting through the neutral 
gas. Under these conditions it may be appropriate to use 
two fluid models which represent the plasma as an ion gas 
interacting with a neutral component (Draine 1980; Toth 
1994; Smith & Mac Low 1997; Stone 1997). 

Recently, Pandey & Wardle (2006) have asserted that 
in weakly ionized plasmas, collisional coupling with the neu- 
trals reduces the effective gyrofrequency of ions by a fac- 
tor pi/ p, where pi is the ion gas density and p is the bulk 
plasma density. The Hall effect then becomes significant un- 
der the relaxed condition L<pc/piU p i. Additionally, given 
the potential importance of charge-carrying grain species in 
molecular clouds (e.g. Wardle 1998, 2004; Ciolek & Roberge 
2002; Falle 2003, hereafter F03), it is clear that a genuinely 
multifluid approach may often be necessary to capture the 
complex interplay of resistive effects due to relative motions 
between species. Similarly, the conditions in proto-stellar ac- 



cretion discs may warrant a multifluid treatment (e.g. Sano 
& Stone 2002a,b; Wardle 2004; Salmeron & Wardle 2005). 

The numerical difficulties introduced by the presence of 
significant Hall diffusion have been outlined by F03, and 
O'Sullivan & Downes (2006, hereafter Paper I). Both of 
these works put forward one-dimensional numerical meth- 
ods for multifluid MHD of weakly ionized plasmas which 
overcome these difficulties. However, the method presented 
in Paper I has the significant advantage of being explicit and 
hence being comparatively easy to implement, particularly 
in codes employing techniques crucial to large scale simula- 
tions such as parallel domain decomposition and adaptive 
mesh refinement (AMR). 

In this paper we present the extension of the method 
described in Paper I to three dimensions. Section 2 details 
the multifluid equations governing weakly ionized plasmas. 
In section 3 we discuss the numerical method used to inte- 
grate these equations, focusing section 3.1 on the treatment 
of magnetic diffusion with particular emphasis on Hall dif- 
fusion. In section 4 we present three-dimensional results of 
shock-tube tests and simulations of three-dimensional turbu- 
lence in both ambipolar and Hall diffusion regimes. Finally, 
in section 5 we make some concluding remarks. 



2 THE MULTIFLUID EQUATIONS 

We assume a weakly ionized plasma such that the mass 
density is dominated by the neutral component of the gas. 
Then, relative to the scale length of the system, if particles 
of a given charged species have small mean free paths in the 
neutral gas, or small Larmor radii, their pressure and inertia 
may be neglected. See F03 for a more detailed discussion. 

For convenience it is assumed that there is no mass 
transfer between species. It is straightforward, however, to 
insert the necessary terms for a more general treatment to 
include mass transfer if necessary. The equations governing 
the evolution of the weakly ionized plasma can then be writ- 
ten as 

%+V-(p n qJ=0, (2) 



-^ + V-(p 1 q 1 q 1 + Pl \)=JxB, (3) 

F) N 

^ + V .[(e 1+Pl )9i] = J-E + J2 H "> W 

n=l 

1 - + V.( ?1 B-B? 1 ) = -VxB', (5) 



a n p n {E + q n x B) + p n p 1 K nl (q 1 - q n ) = 0, (6) 



H n + G„ i + a n p n q n ■ E = 0, (7) 



V • B = 0, (8) 
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J = VxB, (9) 

JV 

^a n p n = Q, (10) 

n=2 
JV 
n=2 

In the preceding equations, the subscripts denote the 
species, with a subscript of 1 indicating the neutral fluid. 
The variables p„, q n = (u n , v n , w n ) T , p n and e„ are the 
mass density, velocity, pressure and total energy respectively 
of species n. In general we assume a closure relation 



E A = -(J x oa) x a A - (18) 
We use the definition 

a x = fxB, (19) 
where X is one of O, H, or A and 

fo = V^E/B, (20) 

/h = ra/B, (21) 

/a = v/^a/- 8 - (22) 



Here ro, ru and ta are the Ohmic, Hall and ambipolar re- 
sistivities respectively defined by the relations 



In 



1 2 



(12) 



where y n is the ratio of specific heats for species n. How- 
ever, for the test cases described here, an isothermal equa- 
tion of state is assumed allowing us to disregard equa- 
tions (4) and (7) and use the closure relation 



: Pi/Pi, 



(13) 



where a is the (constant) isothermal soundspeed. The iden- 
tity matrix, current density and magnetic flux density are 
represented by I, J, B respectively. E' is related to the full 
electric field E by 



E = —q 1 x B + E'. 



(14) 



Additionally, with reference to species n: K n i describes 
the collisional interaction with the neutral fluid, a„ is the 
charge-to-mass ratio, G„ i is the energy transfer rate to the 
neutral fluid, and H n is the energy source or sink. Note 
that in general K n i and G„ i may depend on the tem- 
peratures and relative velocities of the interacting species. 
Equations (2) to (7) are derived from the conservation equa- 
tions for mass (of all species), neutral species momentum, 
neutral species energy, magnetic flux, charged species mo- 
mentum, and charged species energy respectively. Equa- 
tions (8) to (11) describe the solenoidal condition, Ampere's 
law (with displacement current neglected), charge neutral- 
ity, and charge current respectively. We refer the reader to 
F03 and Ciolek & Roberge (2002) for a more detailed dis- 
cussion. 

For a weakly ionized plasma, the generalized Ohm's law 
can be written in terms of contributions from Ohmic, Hall, 
and ambipolar terms (e.g. Wardle & Ng 1999) as 



where 



E' = Eq + En + Ea, 



Eq = (J • ao)ao, 



Ei 



J X ClH, 



(15) 



(16) 



(17) 



ro = 
= 

TA = 



1 



9 i 9 ' 



OA 



with the conductivities given by 



(TO 



1 

- ^2 a n p n f3 n , 



B 



71 = 2 

N 



-y- 

B 1 



a„p„ 



OA 



J_ \ - a n p n 
1 + , 

n = 2 ' 

The Hall parameter /3 n for species n is 

a n B 



Pn 



KlnPl 



(23) 
(24) 
(25) 



(26) 
(27) 
(28) 

(29) 



3 NUMERICAL METHOD 

We assume a piecewise constant solution on a uniform 
mesh of spacing h in each of the x, y and z directions. 
If the solution has been marched forward in time through 
I (not necessarily uniform) intervals, we denote the cur- 
rent time as t and seek the solution at some later time 
t l+1 = t l + t. Cell (i, j, k) of the mesh is defined as the vol- 
ume {(x,y,z) : (i-l/2)/l < x < (i + l/2)h, (j-l/2)h < y < 
(j + l/2)h, (k - l/2)h < z < (k + l/2)h}. Then given any 
quantity D(x, y, z, t) continuously defined on the mesh vol- 
ume, the average value over the cell (i,j,k) at time t l , is de- 
noted by D | ■ k and defined at the cell centre. Note that for 
the sake of clarity we may drop any of the indices i, j , k, or I 
if no ambiguity arises. 

To obtain the full solution at time t l+1 , standard finite 
volume integration methods are applied to all terms in the 
partial differential equations (2) to (5) with the exception 
of the diffusive term —V x E' on the right hand side of 
equation (5) which we discuss in the next section. The time 
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integration is multiplicatively operator split with each oper- 
ation carried out to second order spatial and temporal accu- 
racy in a straightforward extension of the methods described 
in Paper I. Overall second order accuracy in time is main- 
tained by permuting the order of operations (Strang 1968). 
Charged species velocities and pressures may be derived al- 
gebraically by means of equations (6) and (7); the approach 
to the charged velocity is described in appendix A. Finally, 
the V • B = constraint is applied during each timestep. 
Further discussion of this is deferred to section 3.2. 



3.1 Treatment of magnetic diffusion 

We now focus on the numerical methods for integration of 
the magnetic diffusion terms. The induction equation with- 
out the hyperbolic terms is 

= -Vx(£o+£h + £a), (30) 
using equation (15). To proceed, we carry out the expansions 



V x E x 



Fx + F x , 



(31) 



where the subscript X is one of O, H, or A. The correspond- 
ing linear and second order terms, F x and Fx respectively, 
are 



Fo = -[no • (V x J)]oo + [(oo • V) J)] x a 
+a Q V x J, 



(32) 



-[oo ■ (V x a )]J + [(J ■ V)ao] x oo 
+2(J-ao)[Vxao], (33) 



F\ are small in comparison with F , F h , and F\ respec- 
tively. Additionally, under the often reasonable assumption 
that collisional drag on charged particles is dominated by 
magnetic forces, the Ohmic resistivity ro is weak (F03) and 
hence Fq is also small. The stability of a scheme can then 
be investigated through analysis of the reduced induction 
equation 



dB ~ F 1 + F 1 
— — « t H + r A . 

at 



(38) 



The relative importance of the ambipolar and Hall re- 
sistivities may now be parametrized by 77 = rA/|fH|- From 
this point, time intervals are normalized such that f = t/t ± , 
where r x is the characteristic cell crossing time for diffusion 
perpendicular to the magnetic field given by 



h 2 



Equation (38) can be rewritten as 



dB 

dt 



GB, 



(39) 



(40) 



where, using b = B/B, the matrix operator G is given by 
G = Gh + Ga with 



G H = -rn(p- V)(V x •), 

G A = r A [b • (V X (V x •))]& 
-r A [{b- V)(V x ■)] x b. 



(41) 



(42) 



The discretized form of the operator G at time level I, 
denoted G l , is obtained by using the second order derivative 
dicretizations 



F H = (oh ■ V) J, (34) 

f h — -(J-V)o h + (V-oh)J, (35) 

F A = [oa ■ (V x J)]oa - [(o A • V) J)] x oa, (36) 

F A = +[a A ■ (V x a A )]J- [(J- V)a A )] x a A 

-2(J • oa)[V x o A ] + (Va|) x J. (37) 

In the following, we treat the discretization of equa- 
tion (30) as a two part process. Firstly, under certain as- 
sumed conditions, the stability properties of schemes for the 
dominant terms are explored. Secondly, a correction must 
be made to the field updated through such a scheme to in- 
clude any neglected small terms. The latter step is essential 
for consistency with the governing equation and is discussed 
in more detail in section 3.1.3. However for now, we focus 
on the first step of the process. 

Under the assumption of small perturbations in B 
about a mean field, the second order terms F Q , F H , and 



\dx 2 ) t ~ h 2 ' ( ' 

( d 2 B \ _ Bi+ij+i — -Bi+ij-i — Bi-ij+i + Bi-ij-i 
Kdxdy)^ ~ Ih 2 ' 

(44) 

and similar expressions for other terms. Note that schemes 
with simpler discretizations and superior formal stability 
properties may be derived by replacing equation (43) with 
(d 2 B/dx 2 ) t = (B i+2 - 2B t - B,- 2 )/4:h 2 . We do not con- 
sider such schemes further as they are odd-even decoupled 
and hence subject to instability. 

For the purpose for stability analysis, we take a numer- 
ical wave of the form 

B\ jh =B i"-\ (45) 

where Bo is the wave amplitude, i = y/— 1, i = (i, j, k), and 
u> = (uj x , u>y, u> z ). Second order derivatives of B may now 
be replaced using 

d 2 

— — — > \ xx = -2(1 - cos u) x ), (46) 



Modelling weakly ionized plasmas 5 



Hall diffusion 



dx dy 



sm w, sm GJu 



(47) 



and similar substitutions for other terms. A matrix A can 
then be defined whose {x, y) member is given by \ xy . 

Applying the substitutions given by equa- 
tions (46) and (47) to the discretized operators G H and C l A 
yields the skew-symmetric matrix 



(48) 





' 




A H = | 









v Cv" 





and the symmetric matrix 



A A = bC + Cb - tr(A)bb - b 1 £1, 



(49) 



respectively, where £ = l\b, and bC is the dyadic formed from 
b and £. 

With these representations in place, we now look at 
the stability properties of various discretization schemes. 



3.1.1 Standard discretization 

The standard discretization scheme can be written as 



fl+i 



= (I - tG h - tG 1 a )B 1 . 



(50) 



Inserting the numerical wave of equation (45) then 
yields 



B l+1 = (I - ar H A H - ar A A A )B l 

where a = r/h 2 . 



(51) 



Ambipolar diffusion 

Neglecting Ah from equation (51), the eigenvalues of the 
evolution operator (I — ot-aAa) are 



fJ,i = 1 + ar A b £, 
Ma, 3 = l + iar A [tr(A)± |tr(A)6-2C|]. 



(52) 
(53) 



Considering ambipolar diffusion alone, a maximum 
value in the eigenvalue magnitudes is found at a; = 7r(l, 1, 1) 
for an arbitrary orientation of B. The resulting stability 
limit is 



-std , i yj + v 2 

2 77 



(54) 



which is half the corresponding limit for the one-dimensional 
case (Paper I). 



Now neglecting Aa from equation (51), the evolution oper- 
ator (I — qj'hAh) has eigenvalues 



Mi = 

M2,3 = 



1 ± iothC- 



(55) 
(56) 



Clearly, \fj,2, a\ > 1 for all r > 0. The scheme therefore 
requires a vanishing timestep as the Hall resistivity becomes 
large with respect to the ambipolar resistivity such that, as 
in the one-dimensional case (Paper I), 



-STD 
T"H 



as 77 — > 0. 



(57) 



The standard discretization is therefore impractical for 
systems in which the Hall term is dominant. 



Mixed Diffusion 

Equation (51) does not readily allow derivation of general 
analytic expressions for the eigenvalues of the full amplifica- 
tion matrix. However, from the preceding discussions of the 
limiting cases where Hall and ambipolar diffusion terms are 
alternately neglected, and from numerical investigations of 
the intermediate regime, we infer a general case maximum 
in the magnitudes of the eigenvalues when b = (1, 1, l)/v3 
and uj = w(l, 1, 1). Under these assumptions the general 
eigenvalues of the system are 



/Ui = 1 — 2ota(1 — coso;) , (58) 
^12,3= 1 — 2a(rA =F w*h)(1 — costj)(2 + cosa;). (59) 

As rj becomes small, the stability limit is dictated by 
/i2, 3 with a maximum at tu = 2n/3. The corresponding 
timestep limit 



-STD ^ 8 

r < -- 



9 yr 



(60) 



is slightly below the one dimensional limit rj/y/l + rj 2 (Pa- 
per I) and goes to zero with 77. Again, we conclude that the 
standard discretization is impractical for systems in which 
the Hall effect is large. 



3.1.2 Super TimeStepping/Hall Diffusion Scheme 

We now present a technique for overcoming the weaknesses 
of the standard discretization. Similarly to the strategy 
described in Paper I, the induction equation is integrated 
in two parts by multiplicativcly operator splitting the 
Hall and ambipolar terms. A technique known as Super 
TimeStepping (STS) is used to accelerate the timestepping 
for the standard discretization with ambipolar resistivity 
alone. However, STS does not perform well for evolution 
operators with complex eigenvalues and it is evident from 
equation (59) that, for non-zero th and some orientations 
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of b, the eigenvalues may be complex 2 . The Hall term is 
applied separately using a three-dimensional extension of 
the Hall Diffusion Scheme (HDS) introduced in Paper I. 



Super TimeStepping 

STS is a technique which can be used to accelerate explicit 
schemes for parabolic problems. Essentially a Runge-Kutta- 
Chebyshev method, it has been known for some time (Alex- 
iades, Amiez & Gremaud 1996), although it remains poorly 
known in computational astrophysics. 

In this method a "superstep", r STS , is a composite 
timestep built up from a series of Asts substeps such that 



jV STS 



(61) 



Optimal values for dr., yield stability for the superstep 
while the normal stability restrictions on the individual sub- 
steps are relaxed (Alexiades, Amiez & Gremaud 1996). In- 
tegrating the ambipolar diffusion term in this way yields a 
stability limit 



not find any evidence of this in the tests carried out here. 
Under certain conditions however, such as when there is a 
strong directional bias in the initial state, permutation of the 
order may be necessary over successive steps. We anticipate 
such permutation to result in a small reduction in stability 
however. As evidence of this, in the one-dimensional case 
described in Paper I, it can easily be shown that the stable 
timestep limit decreases by a factor of two when the order 
of component updates is alternated. 

In matrix form we can write three-dimensional HDS as 



B l+1 = (l-arHkkA H )(l-arHjjAH)(l-arHnA H )B ! , (67) 

where ii, jj, and kk are dyadics formed from the unit vectors 
i, j, k in the x, y, z coordinate directions respectively. Then 
the eigenvectors of the evolution operator on the right hand 
side of equation (67) are 



Mi = 

M2,3 = 

where 



1- ^±2^9(9-4), 



(68) 
(69) 



-STS -STD N (1 + *Jv) 



(i - V^) 2 



207 (1 + y/v)™ + (1 - 07) 2 



(62) 



(temporarily dropping the STS subscript from N for clarity) 
where v is a user-tunable damping factor and 



,. -STS » r 2 -STD 

hmr A -> A sts t a ■ 

f—*0 



(63) 



STS is first order accurate in time. In order to achieve 
second order accuracy Richardson extrapolation is used. 



g = (ar H ) 2 (C 2 - arnCxty(z). 



Hence for stability we require 



< g < 4. 



(70) 



(71) 



The most stringent restriction is obtained from 
b = (1/a/3)(1, 1, 1) with w = (2tt/3)(1, 1, 1) and related 
symmetry points. Making the appropriate substitutions, 
and additionally using ordinary (unaccelerated) substepping 
with ./Vhds substeps per full timestep, we find 



Hall Diffusion Scheme 

Gjj is skew-symmetric and hence, dropping the H subscript 
for clarity, we can write three-dimensional HDS as 

B l x +1 = B l x -T(G l xv B l y + G l mz B l z ), (64) 
B l y +1 = By — t{G 1 vz B\ + G yx B x +1 ), (65) 
B[ +1 = Bl-r(G l zx B l x +1 + G l zy B l + 1 ). (66) 

Note that equations (64) to (66) are strictly explicit, assum- 
ing they are applied in the order shown, in the sense that 
all terms on the right hand sides are known. However, both 
equations (65) and (66) have implicit- like terms at time t l+1 
on their right hand sides. These terms are the origin of the 
superior stability properties of HDS. 

The order for updating the magnetic field components 
in equations (64) to (66) has been arbitrarily selected. While 
this introduces a directional bias into the scheme, we do 

2 In the one-dimensional case outlined in Paper I, the orienta- 
tion of the field was taken into account explicitly. This allowed a 
finite Hall diffusion term to be admitted while maintaining real 
eigenvalues. 



r^^iVHDS^yiT^, (72) 
V 27 

which is 4/\/27 times the equivalent one-dimensional 
limit (Paper I). Similarly to STS, Richardson extrapolation 
is required to bring HDS to second order temporal accuracy. 

Stability of STS/HDS 

The effective stable timestep limit for the integration of both 
diffusion terms using STS/HDS methods may be estimated 
as the minimum of f™ s and ta TS 

-Sts/hds _ f Th DS if rj <= rf - . 

\ ff TS otherwise, 

where 77* is the solution of f^ 135 — rf TS and depends on the 
user-defined parameters v, Nsts and -/Vhds- 

In the special limiting case given by v = 0, Asts = 1 
and -/Vhds = 1, we have rf = y/27/8. Figure 1 illustrates 
that the stable timestep limit f has a maximum value of 
•\/91/108 at 77 = rj* in this case. The contrast between the 
maximum and minimum possible values of f is then only 
\/91/27. Importantly, f converges to 4/V27 as 77 approaches 
zero unlike the standard scheme for which f goes to zero. 
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Figure 1. STS/HDS for v = 0, AT STS = 1 and 7V H ds = 1: The 
stable timestep limits for HDS (th; solid line) and STS (ta; 
dashed line) as functions of n = r\/\rn\. 



3.1.3 Correction terms 

In the preceding sections, we considered schemes for the ap- 
proximate induction equation (38) in the limit of small per- 
turbations of B about a mean field and small Ohmic resistiv- 
ity ro- As previously stated, however, for consistency with 
equation (30), the neglected small terms must be included 
in the scheme during each update by making the correction 



,1+1 



,z+i 



-riF^+Fl+Fl + F] 



(74) 



All terms are evaluated according to the prescriptions given 
by equations (43) and (44) with the charge current J eval- 
uated via equation (9). 



3.2 V-B = 

It is well known by now that the solenoidal condition on 
the magnetic field is a sensitive issue in any MHD code. 
In our case, however, we have found it to be particularly 
problematic for the tests considered here. 

Both the often inaccurate Powell method (Powell 1994; 
Toth 2000) and the superior Dedner method (Dedner et al. 
2002) rely on reducing the influence of numerically gener- 
ated monopoles by advecting them out of the system, and 
also dissipating them in the case of the Dedner approach. In 
both cases we find that the error, while not fatal, prevents 
convergence in the solution at the expected rate in shock- 
tube tests. Additionally, when periodic boundary conditions 
are employed, advection cannot remove monopoles from the 
system and only the dissipation mechanism of the Dedner 
method has significant effect. 

We find a variant of the constrained-transport (CT) 
method (Evans & Hawley 1988), as described in sec- 
tion 3.2.1, to be effective for periodic boundary conditions 
but impossible to implement with fixed boundary conditions 
in such a way as to obtain a convergent solution for shock- 
tube tests. Fortunately, it is trivial to implement a projec- 
tion technique in this special case as we shall discuss in sec- 
tion 3.2.2. 



3.2.1 Field-interpolated centred differencing 

The family of CT schemes maintain V ■ B by using the 
induction equation to correct the magnetic field generated 
by some base scheme. Usually, this has been done by con- 
structing the electric field on a staggered mesh centred on 
the cell edges. Toth (2000) demonstrates, however, that the 
staggered mesh is unnecessary if a centred differencing of 
the induction equation is carried out on the original grid. 
We make use of the field-interpolated centred differencing 
(field-CD) scheme he presents which has the advantage of 
not requiring any spatial interpolation. 

Field-CD operates by evaluating the electric field E 
on cell centres from the base scheme using the generalized 
Ohm's law given by equation (14). The corrected magnetic 
field B is then given by a centred differencing of the induc- 
tion equation 



B 



i+l 

x i j k 



B 



l 

x i j k 



E 



2h 



^zij-lk 



yij fc+1 



i + l 
yijh- 



(75) 



and similar expressions for the remaining components of B. 

In our case, since we update the magnetic field in an 
operator split fashion, a field-CD correction is made as each 
component of the electric field is applied through the base 
scheme. We find this is more stable than making a single 
correction at the end of a full update via the base scheme. 

Assuming the field is initially divergence- free, equa- 
tion (75) will conserve a centred difference definition of the 
magnetic field divergence 



(V • B)ij k 

Byjj + lk — Byjj-lk 

2h 



+ 



B x j-j-l j J3 X i — lj k 

2h 

2h ' 



(76) 



as long as boundary conditions are compatible. Fixed 
boundary conditions, as required by shock tube tests, are 
not compatible however and an alternative approach must 
be taken. 



3.2.2 Projection 

Projection (Brackbill & Barnes 1980), similarly to CT meth- 
ods, relies on a correction to the magnetic field generated by 
a base scheme. Briefly, the non-solenoidal component of B 
is projected out of the field by solving 



for (f> and making the correction 

B -> B - V0. 



(77) 



(78) 



In Fourier space, writing B = ~^2 m B m , this amounts to 

projecting out the component of each mode B m = e 1 ^ Wm ' r ' 
parallel to the corresponding wavevector u} m using 



B, 



(W • B r 



(79) 
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Figure 2. Neutral fluid x-velocity and {/-component of magnetic 
field for case A with h = 5 X 10 — 3 . The solution from the steady 
state equations, as a line, is overplotted with points from the 
dynamic code. 



4 TESTS 

Similarly to F03 and Paper I, we test the numerical al- 
gorithms outlined here against the multifiuid equations 
for weakly ionized gases in the isothermal limit with two 
charged species. 

4.1 Shock tube tests 

Using analytical solutions to one-dimensional problems for 
comparison, we run the tests obliquely to the coordinate 
axes in the (1, 1, 1) direction. An A 3 grid is allocated for 
each problem but the solution is only calculated in a narrow 
beam with a radius of one cell and a finite length such that 
it is contained completely within the grid. All cells external 
to the beam are referenced by their parallel displacement 
along the beam and treated as boundary cells. For parallel 
displacements outside the range of the beam, the cells are 
set to fixed values. Inside the beam a single reference cell is 
chosen at each unique value of displacement and all external 
cells with the same value are duplicated from this cell. In this 
way a properly three-dimensional problem is possible with 
computation only required on a small fraction of the full A 3 
domain. 




0.2 0.4 0.6 0.8 1 1.2 

X 



Figure 3. Neutral fluid re-velocity and {/-component of magnetic 
field for case B with h = 2 X 10 — 3 . The solution from the steady 
state equations, as a line, is overplotted with points from the 
dynamic code. 



Since for this case we know that w = (1, 1, I)/ \/3, equa- 
tion (79) simply says that for the projection method of diver- 
gence cleaning, the longitudinal component of the magnetic 
field must be held constant as expected trivially from the 
one-dimensional analogs of the solenoidal condition (8) and 
induction equation (5). 

Similarly to F03 and Paper I, the dynamic algorithm de- 
scribed here is tested against solutions of the steady isother- 
mal multifiuid equations. These steady state equations are 
solved using an independent code. The conditions for each 
of the tests are given in Table 1, including the user-defined 
parameters v, Asts and Ahds for STS/HDS substepping. 



4-1.1 Case A: Ambipolar Dominated 

In this test r = 2 x lCT 12 , r H = 1.16 x 1CT 5 and 
ta = 0.068 giving 77 = 5.86 x 10 3 and hence it can be 
expected that ambipolar diffusion will dominate the solu- 
tion. From equation (62) we estimate an overall speed-up 
of about a factor of 2 in comparison with the standard 
explicit approach. Figure 2 shows plots of the x component 
of the neutral velocity, along with B y for both the dynamic 
and steady state solutions. The calculation shown has 
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Table 1. Test calculation parameters. 




Figure 4. Neutral fluid x-velocity and ^/-component of magnetic 



field for case C a 
state equations, 
dynamic code. 



dth h 1 X .1.0 . The solution from the steady 
as a line, is overplotted with points from the 



h — 5 x 10~ 3 . Clearly the agreement between the two 
solutions is extremely good. 

Since the algorithm is designed to be second order 
it is worthwhile measuring the convergence rate of the 
dynamic solution against the solution from the steady 
state solver. The comparison is made using the LI error 



U-2 




Figure 5. Negatively charged fluid x-velocity for case C with 
h = 1 X 10~ 3 . The solution from the steady state equations, as a 
line, is overplotted with points from the dynamic code. 



norm, ei, between a section of the dynamical solution and 
the steady state solution. Working from the downstream 
side, the section xl < x < xr is fixed about the point x* 
where the deviation from the downstream state first ex- 
ceeds 1% of the maximum variation in the solution. Using 
Xl = x* - 0.2 and xr = x* + 0.8 yields ei = 1.00 x 10~ 5 for 
h = 5 x 10~ 3 , and ei = 9.41 x 10~ 5 for h = 1 x 10" 2 . This 
gives ei oc h 3 2 , above the second order convergence ex- 
pected. This may be because of cross-term cancellations aris- 
ing from symmetry in the (1, 1, 1) choice for the direction 
of variation in the problem. 



4-1.2 Case B: Hall Dominated 

The Hall term dominates in this test such that the overall 
efficiency of the scheme is governed by HDS. The parameters 
chosen are r = 2 x 10~ 9 , r H = 0.0116, r A = 5.44 x 10~ 4 
with T) — 0.046 <S l 3 . From equations (72) and (60) we esti- 



3 If the Hall diffusion is increased much further, it appears that 
the approximation of negligible charged particle inertia breaks 
down. 
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mate the scheme to be approximately 20 times faster than 
the standard explicit case. Figure 3 shows the results of the 
calculations for the test with h = 2 x 1CF 3 . For standard 
explicit codes the conditions lead to prohibitive restrictions 
on the timestep. However, the use of HDS allows us to main- 
tain a timestep close to the Courant limit imposed by the 
hyperbolic terms throughout the calculations. 

As with case A, the dynamic solution is tested to ensure 
it has the correct second order convergence characteristics. 
Setting x* at the point where the solution deviates from 
the downstream state by 10% and using xl = x* — 0.05 and 
x R = x* + 1.0, we find ei = 5.11 x 10~ 3 for ft = 2 x 10" 3 
and ei = 1.83 x 10~ 2 for ft = 4 X 10 -3 , giving ei tx ft 1 - 8 . 
The deviation from second order in this case is due to some 
post-shock noise in the high resolution run. 



4-1.3 Case C: Neutral subshock 

This test is similar to case A, but with a higher soundspeed 
and upstream fast Mach number. As a result, a subshock de- 
velops in the neutral flow because the interactions between 
the charged particles and the neutrals are not strong enough 
to completely smooth out the strong initial discontinuity in 
the neutral flow. The ability of the algorithm described to 
deal with discontinuities in the solution is therefore tested. 
Similarly to case A, we expect an overall speed-up of about 
a factor of 2 in comparison with the standard explicit ap- 
proach. 

Figure 4 shows the results of the calculations for 
ft = 1 x 10~ 3 . The subshock in the neutral flow is clearly 
visible as a discontinuity in m, while there is no corre- 
sponding discontinuity in B y . Figure 5 contains a plot of 
the x component of the velocity of the negatively charged 
fluid. There is no discontinuity in this variable, but there are 
some oscillations at the point where the discontinuity in the 
neutral flow occurs as already commented on by F03 and 
Paper I. 

It can be expected that, since there is a discontinuity 
in the solution of this test and a MUSCL-type approach is 
used, the rate of convergence of the dynamic solution will 
be close to first order, at least for resolutions high enough to 
discern the subshock in the solution. Setting x* at the point 
where the solution deviates from the downstream state by 
1% and using xl = x* — 0.02 and xr — x* + 0.1. We find 
ei = 6.44 x 10~ 3 for ft = 1 x 10~ 3 and ei = 1.16 x 10~ 2 for 
ft = 2 x 10~ 3 yielding ei tx ft , close to the first order rate 
anticipated. As in Paper I, we suggest the deviation from 
first order is due to a discontinuity in the electric field at 
the subshock causing an error in the charged velocities since 
smoothing the solution with artificial viscosity improves con- 
vergence. 

4.2 Three-dimensional MHD turbulence 

We now examine the influence of Hall and ambipolar diffu- 
sion on weakly ionized plasmas under the influence a uni- 
form magnetic field Bo superimposed with a weak turbulent 
spectrum of plane waves. Wardle & Ng (1999) assert that the 
system will evolve quite differently depending on which form 
of diffusion is dominant with direct consequences for molec- 
ular cloud support, angular momentum transport in accre- 
tion discs, and dynamo efficiency (Wardle 1998, 2004, 1999; 



Salmeron & Wardle 2005; Sano & Stone 2002a,b; Mininni, 
Gomez & Mahajan 2005). 



4-2.1 Initial B -field generation 

A turbulent field may be represented in a straightforward 
way as a sum of M Fourier modes as 



Bi(r) = A ^ 



■r+p m )s 

*ST7 



(80) 



where A, /3, u> and £ are the amplitude, phase, wave vector, 
and polarization vector of each mode respectively. In the 
limit of a continuous derivative, the solenoidal condition re- 
quires £ m • u> m = for all values of m, i.e.the magnetic field 
is always perpendicular to the direction of propagation. 

Taking a unit cube of 100 3 cells as the computational 
domain this sets a limit on the maximum allowable wave- 
length of Xmax = l/v3. Furthermore, to ensure all modes 
are properly resolved initially, we set the minimum wave- 
length Xmin to 20% of \ ma x such that there are more 
than 10 cells resolving each cycle. Logarithmic spacing in 
uj is then assumed such that Aw m /w ro is a constant where 
ALu m = w TO +i — io m . The amplitude A(n) of each mode is 
generated by 



where 



{u m Lc) 1 



(81) 



(82) 



The variance of the turbulent field is a =< B\ > 
through which the turbulence level E is defined by 
E = a 2 /(B 2 + a 2 ). In the studies below we shall consider 
E — 0.01 and take the variance of the total field < B 2 > to 
be unity such that the Alfvenic signal speed with respect to 
the mean magnetic field is also unity. The correlation length 
L c is set to be Xmax and the normalization factor AV m for 
three-dimensional turbulence is given by 



(83) 

Lastly, for a three-dimensional Kolmogorov spectrum, we 
use a spectral index F = 11/3. 



4.2.2 Results 

For a first approximation to the field we use, a Mersenne 
twister algorithm 4 is called to generate values for the phase 
flm, the direction of u> m , and the orientation of £ m for 
M = 1000 modes. However, this field is neither divergence 
free nor periodic and must be modified. 

Firstly, to derive a periodic field, the components of tJ m 



4 See www.math.sci.hiroshima-u.ac.jp/ m-mat/MT/cwhat-is- 
mt.html. 
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Figure 6. Density power spectra for the Hall (solid 
curve; 0.653 < p < 1.459) and ambipolar (dashed curve; 
0.891 < p < 1.143) cases. A Kolmogorov power law (solid 
straight line) is also shown for reference. 




Figure 7. Ambipolar model kinetic enstrophy with isosurfaces 
at fi = 0.12 (Qmax = 0.64 and < Q >= 0.06). 



must be integral multiples of 2-n. To achieve this, the compo- 
nents are collapsed onto the closest lower integral multiple. 
Secondly, since our measure of V • Bi is not continuous but 
discrete, the above field will not appear divergence free ini- 
tially. Relaxing the condition £ m ■ w m = and assuming a 
centered difference approximation to V • Bi on a grid of 
uniform spacing h then yields the constraint 

i m ■ sm(u m h) = (84) 

for a numerically divergence free field. Since sin(a>/i) is 
known explicitly, we take the polarization with respect to 
this quantity in order to construct an appropriate £. For the 




0.00 

Figure 8. Hall model kinetic enstrophy with isosurfaces at 
H = 1.2 (H max = 4.8 and < Q >= 0.5). 



particular set of modes generated for these tests, the above 
treatment results in 115 unique wavevectors. 

Once B 1 has been fully specified, the direction of the 
mean field Bo is determined by taking a weighted average 
of the wavevector directions as follows 

Bo - B . (85) 

Z^ m =l 

In this case we find B = (0.686, 0.608, -0.399). 

Starting from an initially uniform plasma, Fig. 6 shows 
the density power spectra after 5 crossing times. Clearly 
there is far more structure at all scales for the Hall case 
(except for some low power grid-scale noise at high frequen- 
cies). This behaviour should have significant consequences 
for any gravitationally unstable system. 

Figure 7 shows isosurfaces of enstrophy, denned as 
fi = |V x qj\' 2 , for the ambipolar test with isosurfaces at 
fi = 0.12 (fi max = 0.639 and < fi >= 5.75 x 10~ 2 ). In this 
case the flow has developed vortex tubes about the mean 
field direction. Figure 8 shows isosurfaces of enstrophy for 
the Hall test with isosurfaces at fi = 1.2 (fi ma x = 4.80 and 
< fi >= 0.462). In this case the flow is more complicated 
showing blobs of high vorticity throughout the domain and 
a total enstrophy almost an order of magnitude greater than 
the ambipolar analog. 

Given its relevance to the study of dynamo action (e.g. 
Mininni, Gomez & Mahajan 2005), we also analyze the mag- 
netic helicity defined by H = A ■ B, where B — V x A. Fig- 
ure 9 illustrates the power spectra for both tests. The mag- 
netic helicity is greater at all scales for the Hall case (again, 
except for some low power grid-scale noise at high frequen- 
cies). Initially, we have V < H 2 > = 0.0216, however, in the 
ambipolar case, by the end of the simulation the helicity has 
been largely dissipated to V < H 2 > = 0.0056. On the other 
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Figure 9. Helicity power spectra for the Hall (solid 
curve; l#l max = 1-81 x 1CT 3 , < \H\ >= 3.08 X 1(T 4 , 

< H >= 3.35 X 10 — 9 ) and ambipolar (dashed 
curve; \H\ max = 6.86 X 10" 3 , < \H\ > = 7.36 X 10~ 4 , 

< H >= -1.05 x 10~ 5 ) cases. A Kolmogorov power law 
(solid straight line) is also shown for reference. 



hand, for the Hall regime test V < H 2 > = 0.0142, show- 
ing helicity is well preserved. Clearly, ambipolar and Hall 
diffusion have dramatically different influences on magnetic 
helicity. 



5 CONCLUSIONS 

We have presented a three-dimensional numerical method 
for integrating the multifluid equations appropriate to 
weakly ionized plasmas. Crucially, the method does not rely 
on implicit solvers to counter the poor stability properties 
of conventional explicit schemes. The problematic V x E' 
term describing magnetic diffusion is split into symmet- 
ric and skew-symmetric components representing ambipo- 
lar and Hall diffusion respectively (plus higher order terms). 
The symmetric ambipolar diffusion operator is accelerated 
via the Super TimeStepping (STS) method and the skew- 
symmetric Hall diffusion operator is treated by means of 
the new Hall Diffusion Scheme (HDS). A notable advantage 
of STS/HDS over the standard discretization is that in the 
limit of pure Hall diffusion, the stable timestep limit does 
not vanish. 

Tests are presented for the special case of an isother- 
mal three-fluid gas. For oblique shock tube problems the 
algorithm is accurate and converges approximately to sec- 
ond order when the solution is smooth and to first order 
when the solution contains a discontinuity. We also present 
simulations of magnetic turbulence in the ambipolar and 
Hall regimes and find that the evolution of the gas is very 
different in each case. This result may have profound im- 
plications for environments such as dense molecular clouds 
where magnetic turbulence is important in supporting the 
cloud against gravitational collapse as well as facilitating the 
formation of dense cores. 

The local nature of the explicit scheme means it is 
straightforward to extend to a parallelised AMR context. 



This is in contrast to implicit methods for which this exten- 
sion is difficult. 
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APPENDIX A: CHARGED VELOCITIES 

For this work the collisional coefficients K n i are assumed to 
be independent of velocities and temperatures. The follow- 
ing derivation (S.A.E.G. Falle, private communication) is a 
simplified version of the procedure outlined in Paper I. 

Transforming to the frame co-moving with the neutral 
gas, equation (6) can be written as 

q' n X B - Kn q' n = ~& (Al) 

where k„ = piK„i/a„ and E' may be derived from equa- 
tions (15) through (29) and equation (9). 

The solutions for the charged species' velocities are 
given by 

q' n = -A- X E' (A2) 

where 

(—Kn B z -By \ 
— B z -Kn B x . (A3) 
By B x Kn J 

As in Paper I, this procedure must be carried out itera- 
tively if the collisional coefficients K n i are in fact dependent 
on the velocities of the charged species. If also required for 
K n i , equation (7) may be used to derive the temperatures. 

We point out that interpolating the primitive quanti- 
ties to the cell edges before calculating the charged veloc- 
ities achieves smoother results than by calculating the ve- 
locities at the cell centres and subsequently interpolating to 
the edges. 



